################################################################################
################################################################################
# Who Gets Flagged? An Experiment on Censorship and Bias in Social Media Reporting
# Feezell, Conroy, Gomez-Aguinaga, Wagner
# Last updated: 10/14/2022
# Original Author: Jessica Feezell
# Edited Version: John K. Wagner
################################################################################
################################################################################

# Libraries
library(haven)
library(tidyverse)
library(ggeffects)
library(effects)
library(lme4)
library(ggpubr)
library(stargazer)

# Load Data
Knight <- read_sav("Knight Foundation_FX_SPSS_PLUS SCHOLAR DATA_091521.sav")

################################################################################
################################################################################
################################################################################
# Data Cleaning & Variable creation for analysis----

# Creating Count of Reports Variables
# Representative is a rapist!
Knight <- 
  mutate(Knight, 
         representative = case_when(Q65D_3 == 1 ~ 1,
                                    Q65D_4 == 1 ~ 2,
                                    Q65D_5 == 1 ~ 3,
                                    Q65D_6 == 1 ~ 4,
                                    Q65D_7 == 1 ~ 5
         ))
Knight$representative <-
  factor(Knight$representative,
         levels = c(1, 2, 3, 4, 5),
         labels = c("Black male", "Black female", "White male", "White female", "Control"))

# Waterboard the terrorists!
Knight <- 
  mutate(Knight, 
         terrorist = case_when      (Q66D_3 == 1 ~ 1,
                                     Q66D_4 == 1 ~ 2,
                                     Q66D_5 == 1 ~ 3,
                                     Q66D_6 == 1 ~ 4,
                                     Q66D_7 == 1 ~ 5
         ))
Knight$terrorist <-
  factor(Knight$terrorist,
         levels = c(1, 2, 3, 4, 5),
         labels = c("Black male", "Black female", "White male", "White female", "Control"))

# Cops are racist!
Knight <- 
  mutate(Knight, 
         cops = case_when          (Q67D_3 == 1 ~ 1,
                                    Q67D_4 == 1 ~ 2,
                                    Q67D_5 == 1 ~ 3,
                                    Q67D_6 == 1 ~ 4,
                                    Q67D_7 == 1 ~ 5
         ))
Knight$cops <-
  factor(Knight$cops,
         levels = c(1, 2, 3, 4, 5),
         labels = c("Black male", "Black female", "White male", "White female", "Control"))

# Antifa should die!
Knight <- 
  mutate(Knight, 
         antifa = case_when        (Q68D_3 == 1 ~ 1,
                                    Q68D_4 == 1 ~ 2,
                                    Q68D_5 == 1 ~ 3,
                                    Q68D_6 == 1 ~ 4,
                                    Q68D_7 == 1 ~ 5
         ))
Knight$antifa <-
  factor(Knight$antifa,
         levels = c(1, 2, 3, 4, 5),
         labels = c("Black male", "Black female", "White male", "White female", "Control"))

# Condition / Reported Variables

# Representative is a rapist!
Knight <- mutate(Knight,
                 Q1condition = case_when(Q65D_3 >= 0 ~ 1,
                                         Q65D_4 >= 0 ~ 2,
                                         Q65D_5 >= 0 ~ 3,
                                         Q65D_6 >= 0 ~ 4,
                                         Q65D_7 >= 0 ~ 5)  )
Knight$Q1condition <-
  factor(Knight$Q1condition,
         levels = c(1, 2, 3, 4, 5),
         labels = c("Black male", "Black female", "White male", "White female", "Control"))
Knight$Q1condition <- relevel(Knight$Q1condition, ref = 5) # set reference cat = Control


Knight <- mutate(Knight,
                 Q1reported = case_when(Q65D_3 == 0 ~ 0,
                                        Q65D_3 == 1 ~ 1, 
                                        Q65D_4 == 0 ~ 0,
                                        
                                        Q65D_4 == 1 ~ 1,
                                        Q65D_5 == 0 ~ 0,
                                        Q65D_5 == 1 ~ 1,
                                        Q65D_6 == 0 ~ 0,
                                        Q65D_6 == 1 ~ 1,
                                        Q65D_7 == 0 ~ 0,
                                        Q65D_7 == 1 ~ 1))

# Waterboard the terrorists!
Knight <- mutate(Knight,
                 Q2condition = case_when(Q66D_3 >= 0 ~ 1,
                                         Q66D_4 >= 0 ~ 2,
                                         Q66D_5 >= 0 ~ 3,
                                         Q66D_6 >= 0 ~ 4,
                                         Q66D_7 >= 0 ~ 5))
Knight$Q2condition <-
  factor(Knight$Q2condition,
         levels = c(1, 2, 3, 4, 5),
         labels = c("Black male", "Black female", "White male", "White female", "Control"))
Knight$Q2condition <- relevel(Knight$Q2condition, ref = 5) # set reference cat = Control


Knight <- mutate(Knight,
                 Q2reported = case_when(Q66D_3 == 0 ~ 0,
                                        Q66D_3 == 1 ~ 1, 
                                        Q66D_4 == 0 ~ 0,
                                        Q66D_4 == 1 ~ 1,
                                        Q66D_5 == 0 ~ 0,
                                        Q66D_5 == 1 ~ 1,
                                        Q66D_6 == 0 ~ 0,
                                        Q66D_6 == 1 ~ 1,
                                        Q66D_7 == 0 ~ 0,
                                        Q66D_7 == 1 ~ 1))


# Cops are racist!

Knight <- mutate(Knight,
                 Q3condition = case_when(Q67D_3 >= 0 ~ 1,
                                         Q67D_4 >= 0 ~ 2,
                                         Q67D_5 >= 0 ~ 3,
                                         Q67D_6 >= 0 ~ 4,
                                         Q67D_7 >= 0 ~ 5)  )
Knight$Q3condition <-
  factor(Knight$Q3condition,
         levels = c(1, 2, 3, 4, 5),
         labels = c("Black male", "Black female", "White male", "White female", "Control"))
Knight$Q3condition <- relevel(Knight$Q3condition, ref = 5) # set reference cat = Control


Knight <- mutate(Knight,
                 Q3reported = case_when(Q67D_3 == 0 ~ 0,
                                        Q67D_3 == 1 ~ 1, 
                                        Q67D_4 == 0 ~ 0,
                                        Q67D_4 == 1 ~ 1,
                                        Q67D_5 == 0 ~ 0,
                                        Q67D_5 == 1 ~ 1,
                                        Q67D_6 == 0 ~ 0,
                                        Q67D_6 == 1 ~ 1,
                                        Q67D_7 == 0 ~ 0,
                                        Q67D_7 == 1 ~ 1))


# Antifa should die!
Knight <- mutate(Knight,
                 Q4condition = case_when(Q68D_3 >= 0 ~ 1,
                                         Q68D_4 >= 0 ~ 2,
                                         Q68D_5 >= 0 ~ 3,
                                         Q68D_6 >= 0 ~ 4,
                                         Q68D_7 >= 0 ~ 5)  )
Knight$Q4condition <-
  factor(Knight$Q4condition,
         levels = c(1, 2, 3, 4, 5),
         labels = c("Black male", "Black female", "White male", "White female", "Control"))
Knight$Q4condition <- relevel(Knight$Q4condition, ref = 5) # set reference cat = Control


Knight <- mutate(Knight,
                 Q4reported = case_when(Q68D_3 == 0 ~ 0,
                                        Q68D_3 == 1 ~ 1, 
                                        Q68D_4 == 0 ~ 0,
                                        Q68D_4 == 1 ~ 1,
                                        Q68D_5 == 0 ~ 0,
                                        Q68D_5 == 1 ~ 1,
                                        Q68D_6 == 0 ~ 0,
                                        Q68D_6 == 1 ~ 1,
                                        Q68D_7 == 0 ~ 0,
                                        Q68D_7 == 1 ~ 1))


# Aggregate condition vars

# Black Male Condition
Knight <- mutate(Knight,
                 Q1_BM = case_when(Q65D_3 == 0 ~ 0,
                                   Q65D_3 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q2_BM = case_when(Q66D_3 == 0 ~ 0,
                                   Q66D_3 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q3_BM = case_when(Q67D_3 == 0 ~ 0,
                                   Q67D_3 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q4_BM = case_when(Q68D_3 == 0 ~ 0,
                                   Q68D_3 == 1 ~ 1))

Knight$BM = 
  rowMeans(Knight[,c("Q1_BM", "Q2_BM", "Q3_BM", "Q4_BM")], na.rm = TRUE)

# Black Female Condition
Knight <- mutate(Knight,
                 Q1_BF = case_when(Q65D_4 == 0 ~ 0,
                                   Q65D_4 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q2_BF = case_when(Q66D_4 == 0 ~ 0,
                                   Q66D_4 == 1 ~ 1))


Knight <- mutate(Knight,
                 Q3_BF = case_when(Q67D_4 == 0 ~ 0,
                                   Q67D_4 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q4_BF = case_when(Q68D_4 == 0 ~ 0,
                                   Q68D_4 == 1 ~ 1))

Knight$BF = 
  rowMeans(Knight[,c("Q1_BF", "Q2_BF", "Q3_BF", "Q4_BF")], na.rm = TRUE)

# White Male Condition
Knight <- mutate(Knight,
                 Q1_WM = case_when(Q65D_5 == 0 ~ 0,
                                   Q65D_5 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q2_WM = case_when(Q66D_5 == 0 ~ 0,
                                   Q66D_5 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q3_WM = case_when(Q67D_5 == 0 ~ 0,
                                   Q67D_5 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q4_WM = case_when(Q68D_5 == 0 ~ 0,
                                   Q68D_5 == 1 ~ 1))

Knight$WM = 
  rowMeans(Knight[,c("Q1_WM", "Q2_WM", "Q3_WM", "Q4_WM")], na.rm = TRUE)

# White Female Condition
Knight <- mutate(Knight,
                 Q1_WF = case_when(Q65D_6 == 0 ~ 0,
                                   Q65D_6 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q2_WF = case_when(Q66D_6 == 0 ~ 0,
                                   Q66D_6 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q3_WF = case_when(Q67D_6 == 0 ~ 0,
                                   Q67D_6 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q4_WF = case_when(Q68D_6 == 0 ~ 0,
                                   Q68D_6 == 1 ~ 1))

Knight$WF = 
  rowMeans(Knight[,c("Q1_WF", "Q2_WF", "Q3_WF", "Q4_WF")], na.rm = TRUE)

# Control (No ID) Category
Knight <- mutate(Knight,
                 Q1_No_ID = case_when(Q65D_7 == 0 ~ 0,
                                      Q65D_7 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q2_No_ID = case_when(Q66D_7 == 0 ~ 0,
                                      Q66D_7 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q3_No_ID = case_when(Q67D_7 == 0 ~ 0,
                                      Q67D_7 == 1 ~ 1))

Knight <- mutate(Knight,
                 Q4_No_ID = case_when(Q68D_7 == 0 ~ 0,
                                      Q68D_7 == 1 ~ 1))

Knight$No_ID = 
  rowMeans(Knight[,c("Q1_No_ID", "Q2_No_ID", "Q3_No_ID", "Q4_No_ID")], na.rm = TRUE)

# Demographic Controls
Knight <- mutate(Knight,
                 Female = case_when(ppgender == 1 ~ 0,
                                    ppgender == 2 ~ 1))
Knight <- mutate(Knight, 
                 Education = case_when(ppeduc5 ==1 ~ 1,
                                       ppeduc5 ==2 ~ 2,
                                       ppeduc5 ==3 ~ 3,
                                       ppeduc5 ==4 ~ 4,
                                       ppeduc5 ==5 ~ 5))
Knight <- mutate(Knight,
                 Age = ppage)

Knight <- mutate(Knight,
                 Income = ppinc7)

# Dummy for race (White = 1, All else = 0)
Knight <- mutate(Knight,
                 White = case_when(ppethm == 1 ~ 1,
                                   ppethm == 2 ~ 0,
                                   ppethm == 3 ~ 0,
                                   ppethm == 4 ~ 0,
                                   ppethm == 5 ~ 0))
Knight$White <-
  factor(Knight$White, levels = c(0,1),
         labels = c("Not White",
                    "White"))
# White / Black / Other
Knight <- mutate(Knight,
                 W_B_Other = case_when(ppethm == 1 ~ 1,
                                       ppethm == 2 ~ 2,
                                       ppethm == 3 ~ 3,
                                       ppethm == 4 ~ 3,
                                       ppethm == 5 ~ 3))
Knight$W_B_Other <-
  factor(Knight$W_B_Other, levels = c(1:3),
         labels = c("White",
                    "Black",
                    "Other"))
# Race categorical
Knight <- mutate(Knight,
                 Race = case_when(ppethm == 1 ~ 1,
                                  ppethm == 2 ~ 2,
                                  ppethm == 3 ~ 3,
                                  ppethm == 4 ~ 4,
                                  ppethm == 5 ~ 5))
Knight$Race <-
  factor(Knight$Race, levels = c(1:5),
         labels = c("White non-Hispanic",
                    "Black",
                    "Other non_Hispanic",
                    "Hispanic",
                    "Multi-racial non-Hispanic"))

Knight <- mutate(Knight,
                 PID = case_when(QPID100 == 1 ~ 1,
                                 QPID100 == 2 ~ 2,
                                 QPID100 == 3 ~ 3,
                                 QPID100 == 4 ~ 3))
Knight$PID <- 
  factor(Knight$PID, levels = c(1:3),
         labels = c("Republican",
                    "Democrat",
                    "Independent / Other"))
Knight$PID <- relevel(Knight$PID, ref = 1) # set reference cat = Republican

################################################################################
################################################################################
################################################################################
# T-tests Aggregated Q1-Q4 t-test ----
# BM v control for Q1-Q4 p = .80
t.test(Knight$BM, Knight$No_ID)

# BF v Control for Q1-Q4 p = .52
t.test(Knight$BF, Knight$No_ID)

# WM v COntrol for Q1-Q4 p = .35
t.test(Knight$WM, Knight$No_ID)

# WF v Control for Q1 - Q4 p = .48
t.test(Knight$WF, Knight$No_ID)

################################################################################
################################################################################
################################################################################
# Plot Means for Aggregated treatment conditions
# FIGURE 1 ----
# Create a Separate dataset (means) that we transform to fit our plotting needs
means =
  Knight %>%
  dplyr::select(BM, BF, WM, WF, No_ID) %>% #We're selecting the variables we want to plot.
  stack() %>% # We're stacking them into a single column (it adds a second column keeping track of the variable of origin [e.g., BM, WF, etc.])
  group_by(ind) %>% # Group by the conditions
  dplyr::summarise(mean.ci = list(mean_ci(values))) %>% # Ask for the mean and CI of each condition
  unnest_wider(mean.ci) # Split the answer into three separate variables (because it originally provides a list of the mean, CI min, and CI max)
means$ind = factor(means$ind, labels = c("BM", "BF", "WM", "WF", "Control"))
ggplot(means, aes(ind, y)) + # Plot by group the means
  geom_point() + # The point for each mean
  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = .1) +  # Add bars for the CI
  ylab("Mean Reported") + # Y-axis label
  xlab("Treatment Condition") + # X-axis label
  ylim(.33, .46) # Y-axis range
ggsave("Figure 1.tiff", device = "tiff", width = 5, height = 3, dpi = 1200, units = "in")

################################################################################
################################################################################
################################################################################
# Logistic Regressions for H1 (Appendix D) ---- 

# Post 1
modelQ1 <-glm(Q1reported ~ Q1condition + Female + W_B_Other + Age + Income + Education + PID, 
              family = binomial,
              data = Knight)
summary(modelQ1)

# Post 2
modelQ2 <-glm(Q2reported ~ Q2condition + Female + W_B_Other + Age + Income + Education + PID, 
              family = binomial,
              data = Knight)
summary(modelQ2)

# Post 3
modelQ3 <-glm(Q3reported ~ Q3condition + Female + W_B_Other + Age + Income + Education + PID, 
              family = binomial,
              data = Knight)
summary(modelQ3)

# Post 4
modelQ4 <-glm(Q4reported ~ Q4condition + Female + W_B_Other + Age + Income + Education + PID, 
              family = binomial,
              data = Knight)
summary(modelQ4)


################################################################################
################################################################################
################################################################################
# Stacking Data for H2 ----

###
# create my subset dataset containing the variables I'll need, and add on the id variable:
subset_DF = Knight %>% tibble() %>% dplyr::select(Q1reported, Q2reported, Q3reported, Q4reported,
                                           Q1condition, Q2condition, Q3condition, Q4condition,
                                           Female, White, Age, Income, Education, PID, W_B_Other) %>%
  mutate(id = row_number())

# split that dataframe in two 
subset_DF_reported <- subset_DF %>% dplyr::select(Q1reported, Q2reported, Q3reported, Q4reported,
                                           Female, White, Age, Income, Education, PID, id, W_B_Other)

subset_DF_condition <- subset_DF %>% dplyr::select(Q1condition, Q2condition, Q3condition, Q4condition, id)

# Then we pivot the two dataframes longer. I also mutate a new question variable for each that will help us put them back together.
subset_long_reported=
  subset_DF_reported %>%
  pivot_longer(cols = ends_with("reported"), names_to = c("QReported"), values_to = "Reported") %>%
  mutate(Q = str_sub(QReported, 2, 2)) %>%
  dplyr::select(-QReported)

subset_long_condition =
  subset_DF_condition %>%
  pivot_longer(cols = ends_with("condition"), names_to = c("QCondition"), values_to = "Condition") %>%
  mutate(Q = str_sub(QCondition, 2, 2)) %>%
  dplyr::select(-QCondition)

# Lastly, we merge the two longer dataframes together by id and question number:
subset_long = subset_long_condition %>% full_join(subset_long_reported, by = c("id", "Q"))

# Generate new variable for subset_long where treatment_gender = Female and treatment_race = Black
# Treatment Female
subset_long <- mutate(subset_long,
                      Treat_Female = case_when(Condition == "Black female" ~ 1,
                                               Condition == "White female" ~ 1,
                                               Condition == "Black male" ~ 0,
                                               Condition == "White male" ~ 0,
                                               Condition == "Control" ~ 0))
subset_long$Treat_Female <- 
  factor(subset_long$Treat_Female, levels = c(0:1),
         labels = c("Male / Control",
                    "Female"))
# Treatment Black
subset_long <- mutate(subset_long,
                      Treat_Black = case_when(Condition == "Black female" ~ 1,
                                              Condition == "Black male" ~ 1,
                                              Condition == "White female" ~ 0,
                                              Condition == "White male" ~ 0,
                                              Condition == "Control" ~ 0))
subset_long$Treat_Black <- 
  factor(subset_long$Treat_Black, levels = c(0:1),
         labels = c("White / Control",
                    "Black"))

subset_long <- mutate (subset_long, 
                       Black = case_when (White == "White" ~ 0,
                                          White == "Not White" ~ 1))
# Recodes of Treatment Condition
subset_long = subset_long %>% 
  mutate(WBO_White = ifelse(W_B_Other == "White", 1, 0)) %>% 
  mutate(WBO_Black = ifelse(W_B_Other == "Black", 1, 0)) %>% 
  mutate(WBO_Other = ifelse(W_B_Other == "Other", 1, 0)) %>% 
  mutate(Race = factor(W_B_Other, levels = c("White", "Black", "Other"))) %>% 
  mutate(Condition = factor(Condition)) %>% 
  mutate(Condi_Gender = factor(case_when(Condition == "Black female" ~ "Treat_Female",
                                         Condition == "White female" ~ "Treat_Female",
                                         Condition == "Black male" ~ "Treat_Male",
                                         Condition == "White male" ~ "Treat_Male",
                                         Condition == "Control" ~ "Control"))) %>%
  mutate(Condi_Race = factor(case_when(Condition == "Black female" ~ "Treat_Black",
                                       Condition == "White female" ~ "Treat_White",
                                       Condition == "Black male" ~ "Treat_Black",
                                       Condition == "White male" ~ "Treat_White",
                                       Condition == "Control" ~ "Control"))) %>%
  mutate(Q = factor(Q)) %>% 
  mutate(Female = factor(Female))

# Change Ref. Categories
subset_long$Condi_Gender = relevel(subset_long$Condi_Gender, "Control")
subset_long$Condi_Race = relevel(subset_long$Condi_Race, "Control")
subset_long$Race = relevel(subset_long$Race, "White")


# Random Effects Logit - H2 --- (Appendix F)

# Modeling Treat_Black x Race (Model 1, Appendix F)
rem.raceinteraction = glmer(Reported ~ Treat_Black * Race + Female + Age + Income + Education + PID + (1|id), data = subset_long, family = "binomial", nAGQ = 0)
summary(rem.raceinteraction)
exp(fixef(rem.raceinteraction))

pr.rem.raceinteraction = ggpredict(rem.raceinteraction , terms = c("Treat_Black", "Race"))
plot(pr.rem.raceinteraction)  + ggtitle("") + ylab("Reported") + xlab("Treatment Condition") +  ylim(.15,.35) + scale_color_manual(labels = c("White", "Black", "Other"),
                                                                                                                     values = c("#0072B2", "red", "#009E73")) 

### LEFT HALF OF FIGURE 2
f1 = plot(pr.rem.raceinteraction, colors = "bw")  + ggtitle("") + ylab("Reported") + xlab("Treatment Condition") +  ylim(.15,.35) + theme(legend.box.background = element_rect(colour = "black"))

# Modeling Treat_Female x Female (Model 2, Appendix F)
rem.feminteraction = glmer(Reported ~ Treat_Female * Female + Race + Age + Income + Education + PID + (1|id), data = subset_long, family = "binomial", nAGQ = 0)
summary(rem.feminteraction)
exp(fixef(rem.feminteraction))

pr.rem.feminteraction = ggpredict(rem.feminteraction , terms = c("Treat_Female", "Female"))
plot(pr.rem.feminteraction) + ggtitle("") + ylab("Reported") + xlab("Treatment Condition") +  ylim(.15,.35) + scale_color_manual(labels = c("Male", "Female"),
                                                                                                                    values = c("#0072B2", "red"))

# Full Table, Appendix F
stargazer(rem.raceinteraction, rem.feminteraction, 
          coef = list(exp(fixef(rem.raceinteraction)), exp(fixef(rem.feminteraction))), # Exponentiated Coeff.
          p.auto = F, # Ensure it doesn't exponentiate p-values
          type = "text") # Full Table -



### RIGHT HALF OF FIGURE 2 (These were placed side-by-side in a TIF in an image-editing software)
f2 = plot(pr.rem.feminteraction,  colors = "bw") + ggtitle("") + ylab("Reported") + xlab("Treatment Condition") +  ylim(.15,.35) + theme(legend.box.background = element_rect(colour = "black"))


################################################################################
################################################################################
################################################################################
# FIGURE 2 ----
fig2 = ggarrange(f1, f2, 
                 ncol = 2, nrow = 1, legend = "top")
fig2
ggsave("Figure 2.tiff", device = "tiff", width = 6, height = 3.25, dpi = 1200, units = "in")


################################################################################
################################################################################
################################################################################
# Controls Summary Table (Appendix C) ----
Controls =
  Knight %>%
  dplyr::select(Age, Income, Education, W_B_Other, PID, Female)  %>% 
  mutate(W_B_Other = as.numeric(W_B_Other)) %>%
  mutate(PID = as.numeric(PID)) %>%
  summarise_all(funs(min = min, 
                     median = median, 
                     max = max,
                     mean = mean, 
                     sd = sd))



if(!require(flextable)) # For outputting tables! 
  install.packages("flextable")
library(flextable)

sum.stats = Knight %>%
  dplyr::select(Q1condition, Q2condition, Q3condition, Q4condition, Age, Income, Education, W_B_Other, PID, Female) %>% # Grab variables included in models and conditions (for future listwise deletion)
  na.omit() %>% # Listwise deletion
  dplyr::select(Age, Income, Education, W_B_Other, PID, Female) %>% # Restrict to variables we want to summarize
  mutate(W_B_Other = as.numeric(W_B_Other)) %>% # Convert Race to Numeric (so things like mean, median can be calculated)
  mutate(PID = as.numeric(PID)) %>% # Convert Party ID to numeric as well
  psych::describe() %>% # Create our table of descriptive
  as_tibble(rownames="rowname") # Convert to a data.frame tibble

sum.tab = sum.stats %>% dplyr::select(rowname, n, mean, sd, min, median, max, range) %>% # Select down to summary stats we want
  flextable() # Create Table for Export


sum.tab = set_formatter(sum.tab, mean = function(x) sprintf("%.02f", x),
                        sd = function(x) sprintf("%.03f", x) ) 
sum.tab = set_header_labels(sum.tab, rowname = '', # Labels for stats
                            n = 'N', 
                            mean = 'Mean',
                            sd = "SD", 
                            min = "Min",
                            median = "Median",
                            max = "Max",
                            range = "Range") 

save_as_docx("Summary Stats" = sum.tab, # Export to word
             path = "SummaryTab.docx")


################################################################################
################################################################################
################################################################################
# Treatment Means by Party ID (Appendix E)---- 
pidmeans =
  Knight %>% 
  dplyr::select(BM, BF, WM, WF, No_ID, PID) %>% 
  pivot_longer(cols = c("BM", "BF", "WM", "WF", "No_ID"),
               names_to = "ind",
               values_to = "values") %>% 
  filter(!is.na(values) & !is.na(PID)) %>%
  group_by(PID, ind) %>% 
  dplyr::summarise(mean.ci = list(mean_ci(values))) %>% # Ask for the mean and CI of each condition
  unnest_wider(mean.ci)

ggplot(pidmeans, aes(x = ind, y = y, group = PID)) + # Plot by group the means
  geom_point(aes(shape = PID), position = position_dodge(.5)) + # The point for each mean
  geom_errorbar(aes(ymin = ymin, ymax = ymax), 
                width = .1, 
                position = position_dodge(.5)) +  # Add bars for the CI
  ylab("Mean Reported") + # Y-axis label
  xlab("Treatment Condition") + # X-axis label
  ylim(.2, .55) + # Y-axis range 
  labs(shape = "Party ID")

ggsave("Figure Appendix E.tiff", device = "tiff", width = 5, height = 3, dpi = 1200, units = "in")


